program VMD !  works with df   based on code drivenDf Runge K.f95 which is in frame of source (phase velocity=0)   

  use fft_wrap, only   : fft_init
  use efield
  use evolve

  use dimensions, only : jx, jz, jv, jDirection

  use parameters, only : phi0, kz, vPhase,  Dxxxx, Dvvvv, &            ! real,     phys
                         ny,  iWaveFrame, iFixPump, iModx, iModz, &    ! integer,  phys
                         ModAmp,  Lx, SpeckleWidth, u, &               ! real,     phys
                         nt, dt, iself, iskip                          ! numerical

  implicit none

  !-- major arrays ---------------------


  real       :: dfu(jz, jv, jx, jDirection)     ! distribution function
  real       :: uxMat(jx, jz, jDirection)       ! velocity, x-component
  real       :: uzMat(jx, jz, jDirection)       ! velocity, z-component
  real       :: Ex2D(jx, jz)                    ! electric field, x-component 
  real       :: Ez2D(jx, jz)                    ! electric field, z-component


  !-- suplemental arrays --------------

  real       :: E0z(jx, jz)                     ! background E-field
  real       :: Ezinternal(jx, jz)              ! = Ez2D - E0z, diagnostics
  real       :: Density2Dreal(jx, jz)           ! temporarily array, diagnostics

  real       :: f0vec(jv)                       ! background distribution, gaussian
  complex    :: fvc0T(jv)                       ! FT of f0vec

  real       :: ExFilter(jx)
  complex    :: xMod(jx), uPerpMod(jx)          

  complex    :: E0C(jz)                         ! background E-field at the center
  complex    :: Phi0C(jz)                       ! center potential of background E-field    
  complex    :: phi0T                           ! complex amplitude of background potential
  complex    :: phiTot, phiInt                  ! diagnostics


  !-- variables  ----------------------- 

  real       :: t                               ! time 
  real       :: Courant_in_x,  Courant_in_z     ! advection Courant numbers 
  real       :: OmegaBounce

  integer    :: ix, iz, iv, iDirection          ! counters in r-space  
  integer    :: it, itMin, itMax                ! iteration count


  data       it, itMin / 0, 1 /

 
!-------------------------------------------------------------------------------------


      call fft_init()

      !-- set parameters, coordinate space, and the filter --

      call vmd_init(ExFilter)                                          ! output


      !-- reading restart file or create new initial conditions --
      
      call vmd_ic(f0vec, dfu, uxMat, E0C, phi0C, phi0T, fvc0T, t)      ! output 

      call vmd_print_output(ExFilter)                                  ! input


      !-- creating bunch of "exp" vectors ---

      call evolve_init()
      

      !-- creating two modulation vectors --

      call vmd_init_modulation(xMod, uPerpMod)                         ! output 


      !-- creating initial density, its FFT and modulated E field --

      call efield_init(ExFilter, xMod, E0C,  &                         ! input
                       E0z )                                           ! output


      !-- compute intitial electric field --

      call efield_compute_E(dfu, ExFilter, &                           ! input 
                            Ez2D, Ex2D,  OmegaBounce)                  ! output


      !-- diagnostics --

      call write_csi()


      !-- write initial electric field --

      call write2D_vs_xz("E0z_initial_vs_xz.dat", E0z)
      call write2D_vs_xz("Ex_initial_vs_xz.dat", Ex2D)
      call write2D_vs_xz("Ez_initial_vs_xz.dat", Ez2D)



      !-- write headers for files recording time evolution  --


      write(*,*)
      if(iself==1)then
         write(*,*)'it/nt            t         Phi Internal   OmegaBounce     potEx         potEz'
      else
         write(*,*)'it/nt          kt        real(Epsilon)   Total Energy     Kinetic Energy'
      endif



      open(3,file="vstime.txt")
      if(iself==1)then
         write(3,*) &
      '%t   PhiInternal  PhiToTal -Phase(phi)  KineticEnergy Ex^2energy  Ez^2energy  OmegaBounce Courant_in_z'
      else
         write(3,*) &
      '%t   1/|Epsilon|  Re(Epsilon)  Im(Epsilon)  Re(Csi) Im(Csi) KineticEnergy  Phi0Energy   TotalEnergy   f^2 ntrapped'
      endif
      close(3)


!-------------------------------------------------------------------------------------
! BEGIN EVOLUTION
!-------------------------------------------------------------------------------------


 
do it=itMin,itMIn + nt -1  ! BEGIN EVOLUTION

   call evolve_F_advect( dfu, &                               ! input/output
                         uxMat, f0vec)                        ! input

   call efield_compute_E( dfu, ExFilter, &                    ! input 
                         Ez2D, Ex2D, OmegaBounce)             ! output

   call evolve_F_accel(  dfu, &                               ! input/output
                         Ez2D, fvc0T)                         ! input

   call evolve_U(        uxMat, uzMat, &                      ! input/output
                         dfu, Ex2D, &                         ! input
                         ExFilter, f0vec, uPerpMod)           ! input

   call efield_evolve_E0( E0C, Phi0C, Phi0T,  &               ! input/output
                         xMod,  &                             ! input
                         E0z)                                 ! output

   t=t+dt

      
   Courant_in_z = maxval(abs(uzMat))*dt/dz
   Courant_in_x = maxval(abs(uxMat))*dt/dx

   !-- diagnostics --

   if(  ((it/iskip)*iskip.eq.it) .or. (it.eq.itMIn + nt -1) .or. (it.eq.itMin)  ) then      

      Ezinternal = Ez2D - E0z

      phiTot = phi0T*kz*kz

      call efield_get_phi_internal(phiInt)                         ! output

      call diagnostics( t, float(it-itMin)/nt, Courant_in_z, &     ! input
                        OmegaBounce, SpeckleWidth,  &              ! input
                        dfu, f0vec, &                              ! input
                        Ezinternal, Ez2D, Ex2D, phiTot, phiInt)    ! input

      call efield_write_vs_time(t)                                 ! input


   end if


      
   if(  Courant_in_z > 1.0 ) then
      write(*,*)   'Courant for z advection =', Courant_in_z, '  is > 1, stop evolution'
      go to 314
   endif
            
   if(  Courant_in_x > 1.0 ) then
      write(*,*)   'Courant for x advection =', Courant_in_x, '  is > 1, stop evolution'
      go to 314
   endif

end do      ! END TEMPORAL LOOP
 
314 continue

     close(122)  ! 'Ez_vs_kxtime.dat'

     open(5, file='output', position='append')
     write(5,*) 'Dvvvv should be > E_max / p_max^3 = ', maxval(Ez2D) / pvec(np/2)**3
     close(5) 

     itMax = it
 
     write(*,*) 'end temporal evolution'  


!-------------------------------------------------------------------------------------
! END TEMPORAL LOOP
!-------------------------------------------------------------------------------------


     Ezinternal = Ez2D - E0z


     !-- writing 2D data as function of (x,z) --


     call write2D_vs_xz("E0z_final_vs_xz.dat",   E0z)

     call write2D_vs_xz("Ez_internal_vs_xz.dat", Ezinternal)
     call write2D_vs_xz("Ez_total_vs_xz.dat",    Ez2D)
     call write2D_vs_xz("Ex_vs_xz.dat",          Ex2D)

     call write2D_vs_xz("uxPlus_vs_xz.dat",  uxMat(:,:,1) - u )
     call write2D_vs_xz("uxMinus_vs_xz.dat", uxMat(:,:,2) + u )
     call write2D_vs_xz("uzPlus_vs_xz.dat",  uzMat(:,:,1) + vPhase )
     call write2D_vs_xz("uzMinus_vs_xz.dat", uzMat(:,:,2) + vPhase )


     !-- write total density --

     do ix=1,jx
        do iz=1,jz
           density2Dreal(ix,iz) = sum(dfu(iz,:,ix,:)*dv) 
        end do
     enddo
     call write2D_vs_xz("dnTotal_vs_xz.dat", Density2Dreal)



     !-- write density by direction -- 

     do iDirection=1,2
     
       do ix=1,jx
         do iz=1,jz
            density2Dreal(ix,iz) = sum(dfu(iz,:,ix,iDirection)*dv) 
         end do
       enddo
     
       if(iDirection==1)then
          call write2D_vs_xz("dnPlus_vs_xz.dat",  Density2Dreal)
       else
          call write2D_vs_xz("dnMinus_vs_xz.dat", Density2Dreal)
       endif
    
     enddo 

 
    !-- writing 2D distribution function, f(z,v) ---  

    call write_fdist(dfu, OmegaBounce)                ! input



    !-- writing some integral quantities ---  


    call write_means(dfu, fvc0T)                     ! input


    !-- write some specrtral spectrum quantities --


    call efield_write_spectra()

    call write_umat_spectral(uzMat, uxMat)           ! input



    !-- write restart file --
  
    write(*,*)'phi0=',phi0,'    kz=',kz,' vPhase=',vPhase,'   Lx=', Lx,'     SpeckleWidth=',SpeckleWidth
    write(*,*)'begin restart file write' 

    open(4,file="restart")
    write(4,*) jz,jv,t,ny,phi0,kz,vPhase,ny,Dxxxx,Dvvvv,iWaveFrame, iFixPump, iModx,iModz, ModAmp, Lx, SpeckleWidth
    write(4,*) E0C, phi0C, phi0T      
    write(4,*) uxMat 
    write(4,*) dfu
    close(4)

    write(*,*)'restart file written'; write(*,*)
      

501   format(  6400(e12.6,2x)  )
       
   stop;

end


!-------------------------------------------------------------------------------------
! subroutine deAlias
!-------------------------------------------------------------------------------------

      
subroutine deAlias(jz,et); implicit none;  integer i,jz;  complex et(jz)
if( jz<6 ) return
       do i=jz/2  - jz/6, jz/2 + 2 + jz/6;   et(i)=(0.0,0.0);  end do
return;  end



!-------------------------------------------------------------------------------------
